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Abstract 


The concept of a GPS receiver as a tracking facility and a gradiometer as a 
separate instrument on a low orbiting platform offers a unique tool to map the 
Earth’s gravitational field with unprecedented accuracies. The former technique 
allows determination of the spacecraft’s ephemeris at any epoch to within 3 to 10 
cm, the latter permits the measurement of the tensor of second order derivatives 
of the gravity field to within 0.01 to 0.0001 Eotvos units depending on the type 
of gradiometer. The first part of this report describes a variety of error sources in 
gradiometry where emphasis is placed on the rotational problem pursuing as well 
a static as a dynamic approach. In the second part, an analytical technique is de- 
scribed and applied for an error analysis of gravity field parameters from gradiometer 
and GPS observation types. Results are discussed for various configurations pro- 
posed on Topex/Poseidon, Gravity Probe-B and Aristoteles, indicating that ”GPS 
only” solutions may be computed up to degree and order 35,55 and 85 respectively, 
whereas a combined GPS /gradiometer experiment on Aristoteles may result in an 
acceptable solution up to degree and order 240. 
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Chapter 1 

Introduction 


Before a gravity gradiometer was considered as an instrument on a spacecraft, Wolff 
(1969) suggested a mission where one satellite tracks another [satellite to satellite 
tracking or SST] with the intention of measuring the Earth’s gravitational field. 
The low-low version of this idea has been demonstrated in the ATS 6 / Apollo-Soyuz 
mission cf (VonBun et al.,1980), whereas an actual dedicated low-low Gravity Re- 
search Mission (GRM), as considered in a proposal of the National Aeronautics and 
Space Administration (NASA) cf (Keating et al.,1986), was never realized. The 
high-low version of SST was successfully demonstrated by ATS 6 tracking GEOS-3 
cf (Hajela,1978). A similar technique is used for solving lunar and planetary gravity 
models where velocity perturbations of orbiters are observed on Earth as a Doppler 
shift in the returned radio signals. 

Gradiometry can be considered as a variation of the low-low version of SST 
realized inside one satellite cf (Rummel,1986). Currently there exists a proposal 
within the European Space Agency [ESA] to launch a gravity gradiometer satellite 
called Aristoteles in the time frame of 1996 to 1998. The mission objectives are to 
measure the Earth’s gravity field to within 5 mgal for gravity anomalies and 10 cm for 
geoid heights with a spacial resolution of 100 km. Aristoteles will be placed in a near 
circular sun- synchronous orbit [I = 96.33°] at a height of 200 km. The spacecraft 
will carry a 0.01 E/ n/Hz 2-axis gradiometer [instrument frame perpendicular to 
the satellite’s velocity vector] and a GPS receiver which should allow instantaneous 
estimates of the position to within the sub-decimeter noise level. 

The concept of GPS as a tracking facility on an orbiting platform has also been 
suggested for TOPEX/Poseidon [Ocean/Topography Experiment] and GP-B [Grav- 
ity probe B, a relativistic experiment]. The advantage of GPS is that continuous 
accurate tracking is made possible, allowing the estimation of positions and veloci- 
ties of the spacecraft at each epoch along the orbit. Currently tracking is performed 
mainly by means of laser and Doppler measurements from ground based stations to 
satellites implying that the orbit is covered only up to a few percent with actual 
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measurements. 

The technique for analyzing gravity field errors from GPS position estimates 
and gradiometer observations reported here has been applied in preliminary studies 
of Aristoteles, cf (ESA, 1989) and (Koop et al.,1989). Initially we started with a 
technique for treating the gradiometer problem developed in (Colombo, 1988) and 
included later the GPS part after ESA decided to consider GPS, instead of PRARE, 
on Aristoteles. The GPS part in the error analysis is somehow similar to the problem 
described by Smith et al. (1988) for GP-B. Unfortunately a complete error analysis 
of Aristoteles was not directly possible with the available techniques since it requires 
consideration of 1) a non-polar inclination, 2) a limited bandwidth of the gradiome- 
ter possibly with a colored noise spectrum and 3) the treatment of the GPS and 
gradiometer aspect simultaneously. 

This was the reason to reformulate the original technique in a ’’frequency like 
approach” in which the observation equations sue considered for lumped coefficients 
in the spectral domain. This has been done for both the gradiometer and GPS 
observation equations thereby avoiding explicit analytical expressions of elements 
in the normal matrix. Thus any frequency dependent behavior of an instrument 
may be modeled by means of an a priori covariance function for the noise in the 
observations. 

The organization of Chapters and Appendices in this report is as follows. Chap- 
ter 2 treats some principles of gradiometry in view of Aristoteles, the nominal orbit 
definition, a variety of error sources such as orbital errors, rotational effects includ- 
ing scale, coupling and non-linearity of the gradiometer and the problem of self 
gravitation. Chapter 3 describes briefly the mathematics behind the error analysis, 
expressions for observation equations and the derivation of normal equations, fol- 
lowed by a discussion of the results for various cases. Finally Chapter 4 contains 
conclusions and recommendations for this technique and for gravity field improve- 
ment in general. Two Appendices discuss some specific problems encountered, most 
of them are not directly related to the actual problem. 
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Chapter 2 

Gradiometry 


In this Chapter we will discuss the principles of gradiometry, the choice of orbits for 
proposed missions and some error sources inherent to the concept of gradiometry. 


2.1 Principles 

On the surface of the Earth, the most straightforward method to detect gravitational 
acceleration is to measure the time needed for a proof-mass [p.m.] to fall from a 
certain height, or to observe the period of oscillation of a pendulum with a certain 
length. Both experiments, in some way applied in gravity meters, have been carried 
out for geodetic and geophysical purposes to investigate the gravity field in most 
parts of the world, cf (Vam'cek and Krakiwsky,1984). 

Unfortunately, in an orbiting spacecraft, both the pendulum and the drop test 
fail since the satellite itself is continuously falling resulting in a gravity-free environ- 
ment inside the spacecraft. In this case the only effect that can be observed is the 
remaining non-conservative force primarily caused by atmospheric drag or radiation 
pressure acting on the spacecraft. A successfully applied technique is to correct con- 
tinuously the orbit of a spacecraft by means of small thrusters in such a way that 
a p.m. remains in the center of mass [c.m.]. The resulting orbit is called drag-free 
and obeys the equations of motion: 


#=vy + / (2.1) 

where x represents the acceleration vector [in an inertial coordinate system], V the 
gravitational potential function and / additional conservative forces. 

What would happen if one deployed an accelerometer, consisting of the ”p.m. 
under suspension type”, at some distance from the c.m. Clearly something would 
be observed since the p.m. in the accelerometer would tend to behave as an individ- 
ual orbiting satellite "falling” in another trajectory than the c.m. of the spacecraft. 
However the suspension mechanism of the accelerometer would continously drive 
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the p.m. back to some defined local origin of the accelerometer and as a result 
one would observe some force acting on the p.m. If one ignores rotational effects 
of the entire spacecraft then this driving force consists of the difference in gravita- 
tional acceleration between the c.m. of the spacecraft and the local origin of the 
accelerometer. Using the same technique, differences in accelerations due to gravity 
could be observed at any position in the spacecraft relative to the c.m. or differences 
in acceleration could be observed between two or more arbitrarily placed accelerom- 
eters in a spacecraft. When attached to some frame, e.g. four accelerometers on a 
base plate or eight accelerometers on the corners of a cube, the instrument is called 
a gradiometer. 

The differences in acceleration observed between accelerometers can be trans- 
lated to second order derivatives of the potential F, ignoring effects due to rotation 
which will be discussed later on in this Chapter. Essentially this translation is a 
direct consequence of the equations of motion eq. (2.1). For two accelerometers at 
the points P and Q we find that: 


lo = 


dV 

dZ{ 


*i\o = 


dV 

dz{ 


A Taylor expansion gives: 



dV 

dX{ 


d 2 V 

dxidzj 


Axj + 0 (Ax 3 j ) 
P 


where Axj = Xj\q — Ab a result: 


d^V 

dx{dxj 



+ 0( Ax)) 


In total one could observe a tensor of 9 elements of the second order derivatives of 
V of which 5 components are independent due to symmetry of the tensor and the 
Laplace condition for the gravitational potential, cf (Rummel,1986). 

While the dimension of acceleration is given in m/s 2 , second order derivatives 
are represented in units of 1/s 2 since differences in accelerations are divided by 
meters. It is customary to work with so-called Eotvos units [E] which have the 
dimension of 10 _9 /s 2 . State-of-art gradiometers can operate at room temperatures 
with an accuracy of 0.01 E/\^Hz cf (Benz et al., 1988). Gradiometers cooled at 
super conducting temperatures of a few degrees Kelvin operate with accuracies of 
0.0001 E/y/Hz as is described in (Morgan and Paik,1988). 
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2.2 Orbit selection 

2.2.1 Orbital height 

Although gradiometers of high accuracies can be built with the present state of 
technology, there is a need for circular orbits at very low altitudes [160 to 200 km]. 
A spherical harmonic expansion describing the gravity field shows a natural damping 
behavior containing a term (a e /r) ,+1 for the potential function where o e represents 
the mean equatorial radius and r the radial distance between the instrument and 
the center of mass of the Earth and l the spherical harmonic degree cf (Heiskanen 
and Moritz, 1967). 

Several error analysis studies, more or less similar to the method applied in this 
report, indicate that gravity fields up to degree and orders around 300 to 500 can 
be estimated from gradiometers with accuracies ranging from 10~ 2 to 10 -4 E/y/Hz 
in orbits ranging from 160 to 200 km. [A summary of these studies can be found in 
Appendix C of (Morgan & Paik, 1988)]. In most error analysis studies the threshold 
for recovery in terms of degree and orders is usually determined by comparing the 
estimated error [here r.m.s.] per coefficient per degree to some a priori assumed signal 
behavior of the Earth’s gravity field such as Kaula’s rule of thumb, cf (Kaula, 1966b). 

2.2.2 Orbital decay 

At a height of 200 km an orbit decays about 7 km per day due to atmospheric 
drag which depends on the intensity of solar radiation and the condition of the 
Earth’s magnetic field. This is shown in figure 2.1 where the height of Aristoteles is 
displayed as a function of time over a 1 day period. The underlying simulation, cf 
(Ridgway,1990), involved an integration of the equations of motion using the GEM- 
T2 gravitational model, cf (Marsh et al., 1986,1989), and the Jacchia 71 model using 
Cd=3.0, a cross sectional area of 2.3 m 2 , mass=1240 kg, Kp=2.2 and F10. 7=120 
(average=137) 10 -22 W/m 2 /Hz. These parameters are chosen according to the 
specifications of Aristoteles as given in (ESA, 1989). 

The heights shown in figure 2.1 are with respect to a mean equatorial radius of 
6378137 m and show oscillations of the order of 10 km due to a small eccentricity 
and Cm short periodic effects. The dashed line in this figure is the result of fitting 
a first degree polynomial through the height curve, indicating a slope of -6.8 km 
per day. For the Aristoteles mission it is foreseen to correct the orbit frequently to 
prevent a mean height below 190 km which, according to the simulation described 
above, could occur within 2 days when starting at a mean height of 200 km. 

2.2.3 Sun synchronous orbits 

Benz et al. (1988) explain the need of a sun synchronous orbit at 200 km height 
for Aristoteles. This constrains the inclination of the orbit to 96.33° which can be 
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time [hour] 


Figure 2*1: Height as a function of time for a simulated Aristoteles arc over a 1 day 
period. The dashed line is the result of fitting a first degree polynomial through the 
height curve, indicating a slope of -6.8 km per day. 


confirmed by computing the secular motion of ft due to the flattening term Cjo 
[=-0.00108263] of the gravity field: 


dft 3 nC 20 a e 

~dt = 2(1 - e 2 ) 2 a 2 


( 2 . 2 ) 


where n — y/jlja? as is shown by Kaula (1966a). One finds that 2ir/ft equals to one 
year, i.e. the rate of ft is such that the orbital plane rotates about the z-axis of the 
Earth as fast as the Earth revolves around the sun. Thus it appears for an observer 
in the satellite that the sun is always in the same position with respect to the orbital 
plane. Furthermore, in the case of Aristoteles, ft is chosen such that the sun line is 
perpendicular to the orbital plane, resulting in a so-called dawn-dusk trajectory. 

Despite this geometry, the observer will also notice some yearly variations in the 
position of the sun due to obliqueness of the Earth’s rotational suds with respect to 
the ecliptic. Nevertheless sun synchronous orbits provide an efficient means of power 
production by means of solar arrays and a minimal effect of thermal and mechanical 
noise due to occultation. Figure 4.2 on page 67 in (Morgan & Paik, 1988) clarifies 
the gravity gradiometer orbital lighting geometry in a sun synchronous orbit. 
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A consequence of the sun synchronous orbit of Aristoteles is the loss of polar 
coverage in an area with a diameter of 2 x 6.33° around both poles whereas, from a 
geodetic point of view, 90° inclination is preferable. Most global gradiometer error 
analysis studies, like those of (Colombo, 1988) and (Rapp, 1988), assume a complete 
coverage or I = 90° whereas this is not likely to be the case in an actual gradiometer 
mission. It will be shown that a polar gap of 12.66° in diameter is resulting in a 
poorly posed problem when one aims for a complete gravity field recovery up to say 
degree and order 360. 


2.3 Error sources 

There are various error sources playing a role in gradiometry. The following sub- 
sections will discuss the influence of orbit errors, rotational accelerations, non- 
conservative forces, misalignment and self gravitation on the gradiometer. 

2.3.1 Orbit errors 

As the gradiometer is observing some or more components of the tensor of second 
order derivatives, an error is introduced due to the fact that the orbit, and therefore 
the position of the instrument at a given epoch, can be modeled only up to a certain 
accuracy. One can only assume that the gradiometer performed its measurements 
on some computed orbit whereas in reality tensor components are observed on the 
actual orbit. Orbit errors are mostly caused by errors in force models [such as 
the gravity field, atmospheric drag, radiation pressure and tidal models] which are 
required for the computation of the trajectory of the satellite. In this section we 
will discuss the relation between those forces and the corresponding perturbations of 
the gradiometer satellite. We will not discuss orbit errors due to a limited tracking 
coverage or problems inherent to certain ground based tracking systems as described 
in (Marsh et al.,ibid) since they fall outside the scope of this study. 

Perturbations in near circular orbits due to disturbing [or unmodeled] forces 
acting on the spacecraft are approximated by the Hill equations which are derived 
in e.g. (Schrama, 1989a): 

f u — u ~ 2 tiqv - 3tiqU 

f v = v + 2 n 0 u (2.3) 

f w = tS + nlw 

where u, v and w represent radial, along- and cross-track components of the orbit 
error, n 0 the mean motion of the spacecraft in a circular reference orbit and where 
/ u , f v and f w symbolize disturbing accelerations acting on the satellite. There exist 
exact solutions of these differential equations which are homogeneous, particular 
non-resonant and resonant. 
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The homogeneous solution of the Hill equations is found by solving eqns. (2.3) 
for f u = f v = f w = 0. This solution describes the effect of initial state vector errors 
on the orbit cf Kaplan (1976): 

u(t) = a u cos nt + 6 U sin nt + c u 

v(t) = a v cos nt + b v sin nt + c v + d v t (2*4) 

w(t) = a w cos nt + b w sin nt 

where the constants a u through b w are defined by the initial position and velocity 
errors at a given reference time. 

The particular solutions describe the case where there are forcing functions, here 
chosen as Fourier series, in the problem. The non-resonant particular solution is 
found by solving: 

P u cos art + Q u sinu>t = u - 2 n 0 v - 3tIqU 

P v cosut + Q v sinut = v + 2n 0 u (2.5) 

P w cos a )t + Q w smart = w + n\w 


where P u through Q w symbolize time independent constants. The solution of this 
system of equations becomes: 


u(t) 

v(t) 

w(t) 


wPd 4“ 2n^P v , 

cos u it H -V- s sin u )t 


u '( n o 


a> 2 ) 


w(ng - w 2 ) 


(3ng + w 2 )P„ + 2n 0 uQ u . . (3ng + a > 2 )Q V - 2n 0 u> J P u . 

57^2 rt COS ut + T Sin (Jjt 

u/ 2 (ng - a? 2 ) w 2 (tIq — a; 2 ) 


(ng - w 2 ) 


cos art + 


K - w 2 ) 


smart 


( 2 . 6 ) 


showing that singularity occurs where the denominator becomes zero which is the 
case when uj = 0oru;=±no. These cases require separate, so-called resonant, 
solutions which are described in more detail in (Schrama,1989a). The non-resonant 
particular solution of the Hill equations behaves as a so-called linear system meaning 
that disturbing force functions (the input of the linear system) and perturbations 

07 

in the orbit (the corresponding output) occur at the same frequency of — Hz. 
Characteristic is the damping behavior with respect to u> of the non-resonant solu- 
tions which is caused by the denominators u>(tiq - u> 2 ), u> 2 (n 2 - a; 2 ) and (n§ - u> 2 ) 
in eqns. (2.6). Therefore orbital perturbations occur mostly in the lower frequency 
band between approximately 0 and 3 cycles per revolution as is confirmed by various 
studies such as an orbit error simulation described in (Schrama, 1989a). 

A similar damping behavior with respect to frequency can be expected in the 
gradiometer error signal caused by orbital perturbations. This effect is approximated 
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by linearizing second order derivatives of the potential function V — \ijr whereby 
t — (x 2 + y 2 + z 2 ) 1 / 2 . The required third order derivatives take the following form: 


d 3 V 

dx{dxjdxk 


=,{ 


3i fc dxi _ 15 XiXjX k ji 


dxj 


+ 


r 5 dxk 




For xi = x, X 2 = y and x 3 = z this results in the following Taylor expansion: 


" v 

v XX 

v 

v xy 

V*Z ' 


' V xx 

v 

r xy 

V xz ' 

Vi 

' Az 

0 

Ax 

Vy X 

Vyy 


= 

Vy X 

Vyy 

y* 

0 

Az 

Ay 


V*y 

v zz 


V zx 

Vzy 

Vzz 


Ax 

Ay 

— 2Az 


(x+Ax f y+Ay,*+Az) 


in which one replaces Ax = w , Ay = v and Az = u. At 200 km altitude the 
term 3/i/r 4 equals approximately 6.4 X 10“ 13 indicating that orbit errors of 

the order of 10 m are required to obtain the 0.01 E level whereas errors of 0.1 m 
correspond to 10~ 4 E . 

Figure 2.2 shows the results of a simulation cf Bettadpur (1990) in which tensor 
components are computed as if they occur on a reference and a perturbed trajectory. 
The amplitude density spectrum of the differences between T uu on both trajectories 
demonstrates that most of the orbit error problem in the gravity gradients occurs 
below 4 cpr for a 0.01 E/\/H z instrument and below 25 cpr for a 0.0001 E/y/Rz 
instrument. A simple, but efficient, way of avoiding the orbit error problem is there- 
fore to filter the lower part of the frequency spectrum in the error analysis. Other 
techniques to treat the effect of orbit errors on gravity gradients are cf (Rummel & 
Colombo, 1985): 

• to consider orbit error free combinations such as 2V XX — V 2Z and 2 Vyy — V zz , 


• to introduce initial state vector components and possibly forcing terms (P u 
through Q w in eq. (2.5)) as unknowns in an estimation problem. [The obser- 
vation equations for this problem are obtadned by substitution of (2.4) and 
(2.6) in (2.7)]. 

In this study the former technique, elimination of the lower part of the spectrum, is 
used to avoid any unnecessary complexity in the error analysis. Various references 
on the orbit error problem in gradiometry can be found in (Rummel, 1986), the 
technique of filtering originates from (Colombo, 1988). 


2.3.2 Rotational effects 

Any rotation of the gradiometer causes disturbing rotational accelerations which are 
observed by the instrument. Therefore in the following two subsections we will dis- 
cuss 1) the effect of angular velocities [w] and accelerations [d;] on the accelerometers 
[a static approach] and 2) the behavior of u? and c b in time [a dynamic approach]. 
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Frequency [cpr] 

Figure 2.2: Amplitude density spectrum of the error in T uu caused by orbital per* 
turbations. On the horizontal axis the frequency is represented in terms of cycles 
per revolution with a resolution of 0.7 cpr. On the y-axis the error is presented in 
terms of E. 

Attitude problem, static approach 

According to Rummel (1986) any accelerometer will measure the total acceleration 
R which is: 

R — Rq -f- cv x Rq -f- u i x R + u x (u? x iZ) (2*8) 

where j? 0 and 7?o describe respectively the acceleration and the velocity of the in- 
strument frame and where R equals to the displacement vector relative to the center 
of mass of the spacecraft. Equation (2.8) can be evaluated for differences in acceler- 
ations which are actually observed in the instrument frame whose origin is located 
at the center of mass of the spacecraft. This results in: 


■ A ' 


' A R x ' 

ar 2 

= (r + n + n 2 ) 

AJt 2 

A R 3 


ar 3 
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where Ai? symbolizes the difference of the displacement vector of two accelerometers 
while T symbolizes the tensor of second order derivatives of V : 


Furthermore 


and 


n 2 



v 

r XX 

v 

v xy 

v KZ ' 

r = 

V yx 

Vyy 

Vy> 


. F - 

V,y 

V zz 


0 


U?2 

li = 

U>3 

0 

-Wl 


-U>2 

Wi 

0 

1 

3 

1 

«5 




V1W2 




w 3 


U?2^3 


W iU7 3 
U? 2 W 3 
-W? 




For A = T + Ji + fi 2 , the tensor which is actually observed by the instrument, we 
find the following relations: 

( 2 . 10 ) 


and 


-(A- A 2 ) 


r + n 2 = -(a + a t ). 


(2.11) 


In principle, for a 3-axis gradiometer, one could obtain J1 by an integration of ji 
with respect to time, cf (Rummel,1986): 


m = f 

Jo 


Sldt + fi 0 


( 2 . 12 ) 


This approach helps to estimate fi 2 in eq. (2.11) thereby separating rotational from 
gravitational effects. Unfortunately this technique can not be applied for a 2-axis 
gradiometer as is the case for Aristoteles. Assuming that axis number 1 is radial, 2 
along track and 3 cross track, heading in the same direction as the angular momen- 
tum vector of the orbit, we find that the following components can be observed: 


An = I'll — (w| + U>|) 

A 33 = I^ — (u > 2 -f u> 2 ) (2.13) 

A13 = 2^ 13 F ^ 31 ) = ^ 13 F u,lU ’ 3 

Assuming that u> 3 [nominally the mean motion when the spacecraft is designed to 
fly in an Earth pointing mode] is far larger than o>i or o; 2 and ignoring all terms 
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containing w 2 since they are estimable by means of eq. (2.10) and (2.12) one can 
linearize eqns. (2.13) as: 

AAn = ATu — 2noAu>3 

AA 33 = Ar 33 — 2u;xAu;i (2*14) 

AAJ 3 = ATi 3 + Au ; 3 + rioAu;i 

From eqns. (2.14) one can conclude for a 0.0lE/\/H z instrument that Au> x and Au> 3 
[here symbolizing the yaw and pitch rate] must be known respectively to 10' 8 and 
5 X 10“° rad/s which poses a severe constraint on the restitution of attitude of the 
instrument, cf (ESA, 1989). 

According to this report (ibid) modern gyroscopes obtain an accuracy of 10~ 7 
rad/s inside the measurement bandwidth of 5 x 10“ 3 to 0.125 Hz which is unfor- 
tunately still a factor 10 to 20 too large for an adequate attitude reconstruction. 
Even the inclusion of a star tracker [1 arc second or 4.8 X 10" 6 radians is feasible 
by current space qualified star trackers] in this configuration wouldn’t help since the 
required 10 -8 rad/s couldn’t be obtained in a 4 second integration period which is 
the sampling time proposed for Aristoteles. 

A possible solution, proposed by Matra Espace cf (ESA, 1989), is to predict 
with existing gravity models values of T x3 to within 0.5 E in order to improve the 
estimation of u>i and u; 3 thereby enabling to derive more accurately T xx and r 33 . 
This seems to be possible since the differences of T uw computed by two existing 
high degree and order gravity models seem to be smaller than the required 0.5 E, 
cf (Schrama, 1989b). However in this technical memorandum (ibid) we warned that 
both models share mostly the same observations [mostly gravity anomalies] so that 
the statistics are obscuring the real accuracy of Ti 3 . With this in mind Matra’s 
proposal leads to a vicious circle where one builds a gradiometer to measure a high 
degree and order gravity field which happens to operate only when an a priori model 
of such a field exists. 

The above mentioned problems were a good reason to consider the attitude 
problem in a dynamic approach [considering differential equations] in an attempt 
to describe A as functions of the time caused by disturbing torques. Such an ap- 
proach explains the behavior of the rotational velocity components in the frequency 
domain as is shown in the following section. 

Attitude problem, dynamical approach 

The Newton-Euler equations take the following form: 

#= -uf x H + T (2.15) 

where Tf is equal to the angular momentum vector, u is a polar vector containing 
angular velocities and T is a torque vector. By definition, H = /a;, where / is a 
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tensor containing the moments of inertia of the body under consideration. As a 
result the Newton- Euler rotational equations take the following form: 

Iuj — — uj X /w -f T ( 2 . 16 ) 


These equations are considered in the case where I represents a diagonal matrix 
containing the principle axes of inertia: 


Ji*i = {^2 ~ + T \ 

I 2&2 — {I 3 “ -fi)^i^3 + T 2 (2.17) 

Iz&z = (-fi *** + T 3 


In our case we know that Aristoteles is in an Earth pointing mode and that the 
torques are caused by a combination of gravitational torques, control torques needed 
for attitude control [momentum wheels], and other torques which are mainly due 
to the atmospheric drag acting on the satellite. Some insight can be gained by 
linearizing the Newton-Euler equations [including the gravitational torques] for small 
rotations assuming a nominal rotation about the 0*3 axis [cross-track axis] of the 
spacecraft. This results in the following system of differential equations as is shown 
in (Morgan and Paik,1988): 

I\0\ — (A + I 2 — ^3) n 0^2 + (^2 " ^3) n o^i + 

^ 2^2 — (T 3 — I 2 — A) n 0^1 + 4(/i — /3 )tIq02 + T 2 (2.18) 

hh = 3(/i - h)nl9 3 + T 3 


Here the variables 0 t symbolize small angles [8+ = a;,*], whereas T{ symbolizes torques 
free of gravitational effects. The particular non- resonant solution of (2.18) is ob- 
tained by assuming that: 


Ti 

h 

T 2 

h 

Tz 

h 


Rx cos/3nof 
R 2 sin/3nof 
R$ cos 


(2.19) 


whereby /3 is determining the frequency in terms of cycles per revolution. The 
non-resonant solutions become: 


*(0 

e 2 (t) 

^ 3(0 


-{? + Q2)Ri + 0PiR* 

n 2 o ((0 2 + Qi)(P 1 + Q i ) + P 2 PiP 2 ) 
~(/ 3 2 + QQJta - PP2R1 
n 2 0 ((P 2 + Qi)(/3 2 + Q 2 ) + P 2 PiP 2 ) 


cos flnot 
sin /3 not 


-R 3 

n 2 o (0 2 + Q 3 ) 


cosflnot 


( 2 . 20 ) 
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whereby 


Ii + h ~ h ~ h ~ h 

p, = — h — ’ Q, = ~ir 

Fl ~ — h — ’ 

3(Wj) 

V3 h 

These solutions describe the behavior of attitude errors as a result of disturbing 
non-gravitational torques [divided by moments of inertia] which converted to rota- 
tional velocities [Au?x = 0 \ and Aoj 3 = 0 3 ] and substituted in eqns. (2.14) result in 
error estimates for AA y. Note that eqns.(2.20a-b) become singular when: 

P 2 = -\{PxPi + Qi + Q 2 ) ± \>J{PxP% + Qi + Q 2) 2 - 4 Q 1 Q 2 (2.21) 

and that (2.20c) becomes singular when: 

P = -Q 3 . (2.22) 

This problem is not considered any further, it results in a resonant set of D.E.’s 
which must be treated separately. 

More important is the result that attitude errors, and thereby angular velocities 
and gradiometer signal errors, are decaying at a rate inversely proportional to the 
frequency of the disturbing torque function, see also eq.(2.20). Therefore one could 
expect that attitude errors are confined to the lower frequency band and that the 
problem could be avoided by means of a high pass filtering technique similar to the 
way orbital errors are treated. This is discussed in Appendix A where the results 
of a simulation of attitude errors for a 10 -2 and a 10 -4 E gradiometer satellite are 
shown. Unfortunately the results indicate that the torques caused by atmospheric 
drag are probably too large for even a 10“ 2 E instrument and that the attitude 
must be reconstructed to higher accuracies than presently suggested for Aristoteles. 
Additional studies are needed to determine whether this is feasible with modern 
processing techniques and or technology resulting in an acceptable solution within 
the budgetary constraints. 

2.3.3 Other effects 

So far the effects of orbit errors and rotations have been discussed. There are several 
other effects causing errors in the gradiometer signal such as misalignment between 
axes, scale-errors, coupling and non-linearity of the accelerometers, including effects 
due to self gravitation [fuel sloshing]. 
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Scale, coupling and non-linearity 

According to Touboul et al. (1990) the acceleration measured along the i-axis 7 ; may 
be described by: 

7; = (1 + ei)Ti + €?Yj -f ^T k 

+ + eyTiTj + gTiTk (2.23) 

+ bias + noise 


where: 

• Tj, Tj and Tk are the sum of all external accelerations projected on the i j and 
k-axes, 

• is a bias term for the i-axis, e.g. due to an electronic or mechanical bias in 
the accelerometer, 

• e- and e* are the coupling terms between the i j and k-axis, due to misalignment 
of the sensitive axis of the accelerometer and the actual frame axis of the 
gradiometer and obliqueness of axes, 

• ej*, t 1 / and ej k are non-linearity terms which are mainly due to defects of sym- 
metry of the electrostatic suspension system around the accelerometer proof 
mass. 

A complete treatment of the ’’scale, coupling non-linearity” problem for Aristoteles 
is described by Touboul et al. (ibid). They mention that in-orbit measurements by 
means of a calibration device are needed to obtain €( at « 10~ 5 in a relative sense 
and e? and at « 10“ 6 rad [which allows atmospheric drag accelerations up to 
5 X 1(T 7 m/s 2 ]. 

The principle danger of ’’scaling, coupling and non-linearity” errors in the gra- 
diometer is that non-conservative external forces, in this case dominated by atmo- 
spheric drag, enter directly in the observed signal. The magnitude and spectral 
behavior of the drag fluctuations are derived from the results obtained from the 
missions of Castor, Atmosphere Explorer-C and Dynamics Explorer 2 which are 
also discussed in (Touboul et al.,ibid). They conclude that the velocity vector of the 
spacecraft must be as perpendicular as possible to the gradiometer plane. This may 
require a so-called yaw- steering mode of Aristoteles compensating for cross- track 

winds near the poles, see also (ESA, 1989). 

Self gravitation 

Figure 2.3 is taken from (Morgan & Paik,1988) and shows the expected measurement 
signal in terms of E due to masses varying from 0.01 to 1000 kg in the range of 0.1 
to 100 meter from the gradiometer. It explains that any gradiometer in a spacecraft 
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MASS M(kg) 

Figure 2.3: Gradient sensitivity level, cf (Morgan and Paik,1988). 

is biased to a certain amount by self-gravitation. More serious are masses vibrating 
at frequencies inside the measurement bandwidth of the gradiometer, implying that 
care must be taken in the design of antennas and other appendages so that the eigen- 
frequencies of vibration are outside the measurement bandwidth or alternatively to 
insure that the magnitude of vibration is too small to be noticed by the instrument. 
In this context sloshing of fuel should be avoided in the proximity of the instrument 
requiring a special design of the hydrazine tanks, see also (ESA, 1989). 
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Chapter 3 

Gravity field error analysis 


3.1 Introduction 

The objective for launching a satellite equipped with a gradiometer and a GPS 
receiver is to improve the Earth’s gravity model. The goal of an error analysis is 
to quantify the expected accuracy of recovered potential coefficients given certain 
characteristics of the instruments, the orbit and a priori information about the 
gravitational field. 

The next logical step after an error analysis would be to carry out an actual sim- 
ulation/recovery experiment. In the simulation part of such an experiment, gravity 
gradients including noise and systematic effects are generated by means of existing 
models and known characteristics of the spacecraft and instruments; during recovery 
one attempts to estimate the potential coefficients from the simulated observations. 

Undoubt ly the latter experiment is more convincing for demonstrating the effi- 
ciency of a processing technique. However it is also far more laborious than the error 
analysis technique described here and therefore an expensive method for studying 
the effect of assumptions made in the generation part. [The generation part would 
for instance depend on the availability of a super computer and a good deal of com- 
puting time since it consists of evaluating spherical harmonic expressions of gravity 
gradients up to / = 360 along a reference orbit.] 

8.1.1 Problem definition 

The problem definition assumed here is shortly summarized as follows: 

• A circular orbit is assumed, orbital decay, eccentricity effects and C 20 short 
periodic oscillations are not considered, 

• The inclination of the orbital plane I is fixed during the mission, any choice 
of I is allowed, 
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• In the nominal orbit the elements ft,u; and M are allowed to drift linearily as 
a result of C 2 o secular gravitational effects, 

• It is assumed that the orbit repeat ratio allows the estimation of a gravity field 
up to a given degree and order, 

• Gradiometer observations as well as position estimates from a GPS receiver 
are used as observation types to model the gravity field, 

• Covariance functions of the above mentioned observation types are formulated 
in the spectral domain and are based upon the accuracy of the instrument, 
sampling time and the mission duration. 


3.2 Error analysis 

3.2.1 Lumped coefficient approach 

In the error analysis so-called lumped coefficient expressions of gravity gradients and 
GPS position estimates are used as observation equations. Therefore we discuss the 
subjects: spherical harmonics along a reference orbit, the requirements for a repeat 
orbit and the gradiometer and GPS observation equations in the form of lumped 
coefficient expressions. 


Spherical harmonics along a circular orbit 

The potential function, including gravity gradients as will be shown later on, ex- 
pressed in spherical harmonics up to degree and order L projected on the nominal 
orbit may be written as a Fourier series: 

L L 

T= A km COS 1pkm + BkmSiMpkm (3.1) 

k=—L m = 0 


where Ak m and Bkm are so-called lumped coefficients related to the original potential 
coefficients as, cf (Schrama, 1989a): 



(3.2) 
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Imin = max(|fc|, m) + <5, 


( 0 when A: - max(|Ar|,m) : even 
1 when k - max(|i|,m) : odd 


and 

V’fcm = V’fcm + i’kmt (3.3) 

where it is assumed that Tpkm 15 determined by J 2 secular effects as is discussed in 
(Kaula, 1966a): 

ifrkm - + M) + m(fl - 6) - ku a + md7 e (3.4) 


Non-overlapping lumped coefficients 

A block diagonal system of observation equations is obtained when lumped coef- 
ficients for arbitrary values of k and m occur at a unique frequency Hz. 

This convenient property will be used throughout this study and is described in this 
section. 

The actual frequency in cycles per revolution /3* m [= fpkm/^o] may be written as: 

(3km = k + m^- - k + (3.5) 

iV 0 iV r 

where k G [-£,£], m G [0 ,£] and {N d , N r } G A/" due to the orbit repeat condition. 
The variables N d and N r are respectively the number of nodal days [27r/d> c seconds] 
and the number of revolutions in a repeat period. In order to prevent overlapping 
of lumped coefficients one has to avoid: 

^ ^ where \N$\ < \N d \ A |7V r *| < \N r \ and {N* d ,N r *} € Af 

which results in the following conditions: 

• N d can be an arbitrary integer 

• N r must be a prime number. 

Secondly one has to avoid: 

Pkimi - /3fc a m, where li/ljAmi / m 2 

which is possible for k, m combinations resulting in terms which are 180° out 
of phase so that /3fc im , = -(3k 2 m 2 - According to eq.(3.5) this is the case when: 


resulting in: 


Nd _ _ + kj 

N r mi + m 2 

It follows directly that overlapping of fik m can be avoided by taking N r greater than 
2 L since max(m! + m 2 ) = 2 L. However overlapping of zonal lumped coefficients can 
not be prevented since it will always occur for fiko and fi-ko- 

Both conditions are fulfilled when there are a sufficient number of revolutions 
in a repeat period [N r > 2 L] while N r is chosen as a prime number. For an actual 
gravity mission with the objective to solve for a gravity field up to L = 360 this 
means that at least 44.3 days or 721 revolutions are needed in a repeat period. 

Gradiometer observation equations 

The conventional way of expressing gravity gradients at a point somewhere in the 
A space is discussed by Rummel (1986). In his lecture notes gravity gradients 
are expressed in terms of derivatives of the potential function with respect to r, <f> 
and A assuming that tensor components are evaluated on a polar orbit. Here these 
expressions can not be used since gravity gradients are required along an inclined 
orbital plane. 

The expressions used here to relate gravity gradients to partial derivatives of 
orbital parameters are discussed in Appendix B. By using eqns. (B.6), (B.7), (B.8) 
and (B.9) one can derive the following H and G terms as they are used in equations 
identical to (3.2) for lumped coefficients of gravity gradients: 


TTUU 
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s T uw 

and T vw there exist different 

expressions for the 


For the tensor components T, 

lumped coefficients due to a i „ n 

(B.15). In this case the lumped coefficients are related to a/ m and fii m as: [• = u, v 
or w, see Appendix B] 
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The G and H terms become 
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(3.14) 


GPS observation equations 

The GPS observation equations are derived from the non-resonant particular solu- 
tion of the Hill equations given by eq. (2.6). In the next step w in this equation is 
replaced by /3km n 0 and all partial derivatives are substituted [• = u, v or w]: 

A • (0 = E £ cos V’fcm + B& sinV’fcm (3.15) 

k m 


The observation equations for u and v become: 
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3.2.2 Normal equations 

The system of observation equations: 

j/ = A® + e 


(3.16) 
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(3.18) 
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is formed from a combination of gradiometer and the GPS observation equations 
derived in the previous sections, [namely eq.(3.2) using Hi m k and terms ac- 

cording to eqns. (3.6) through (3.9) and eq.(3.10) through (3.14) for gradiometry 
and (3.15) through (3.20) for GPS.] [The vector y contains the observations in the 
form of lumped coefficients, the vector x contains unknowns for Ci m and S/ m .] The 
A matrix, containing the partial derivatives of the lumped coefficients with respect 
to the unknowns, is block diagonal [one block per order m] provided that the repeat 
ratio of the reference orbit is chosen such that there is a non overlapping lumped co- 
efficient configuration. In a least squares approach [minimizing e T Q yy t] the normal 
equations become: 

i = (A T Q^A)~ 1 A T Q^y (3.22) 

where N = A T Q^A is called the normal matrix and where Q yy is a covariance 
matrix of the observations which is considered to be diagonal. In this case N is also 
block diagonal so that the algorithm for building up and inverting the entire normal 
matrix is a sequential process in which each block is treated individually. 

It is well known that the inverse of the normal matrix equals the covariance ma- 
trix of the estimated parameters as is discussed in e.g. (Schrama, 1989a). As a result, 
the diagonal elements of N~ l are the estimated variances of potential coefficients 
which are derived from the observation equations used to build the normal matrix. 
In the error analysis described here these diagonal elements are used for computing 
the r.m.s. values per degree of potential coefficients, gravity anomalies and geoid 
heights as will be discussed later on. 


A priori observation variance model 


The diagonal Q yy matrix mentioned in the previous section contains the a priori 
variances of the observations [on the main diagonal] which are in our case the lumped 
coefficients [for GPS position as well as gradiometer observation equations]. The 
variance to assign per lumped coefficient depends on the instrument accuracy < 7 /, 
the sampling time At and the total length in time over which the samples are taken 
T [also referred to as the mission duration]. Here it is assumed that 07 represents a 
r.m.s. value of all samples in the set of observations. By propagation of variances 
one can show that: 

/A t\ l ! 2 

= *1 (=r) (3.23) 


<?o 


where <r 0 equals to the individual r.m.s. per sample. If one assumes that there exists 
a uniform flat noise spectrum for the observations and that lumped coefficients are 
obtained from the Fourier transformation of the observation sequence then a 0 equals 
to the r.m.s. per lumped coefficient [due to Parceval’s identity]. 


Unfortunately, the total noise spectrum for the gradiometer is not flat; instead it is 
band limited from 5 X 10 -3 to 0.125 Hz, cf (ESA, 1989), meaning that (3.23) may 
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be applied only inside this frequency band. In order to simulate the attitude error 
problem a 1//3 behavior is assumed for the r.m.s. between 4 and 27 cpr, [/3 m ; n = 27 
cpr corresponds to 5 x 10~ 3 Hz.] Accordingly the r.m.s. below 27 cpr is modeled as 
(Anm//3) x <To; below 4 cpr an infinite r.m.s. is assumed to avoid the ”orbit problem” 
in the error analysis [see Chapter 2], 

In a so-called best case analysis the noise spectrum of GPS position observations is 
considered to be flat using eq.(3.23). In the worst case Smith et al. (1988) mention 
that frequencies which are multiples of once per revolution modulated by multiples 
of once per 12 hours again modulated by multiples of once per day, are omitted in the 
observation noise spectrum. The rationale is that l) frequencies which are multiples 
of once per revolution are caused by failing to solve properly for the trajectory of 
the low orbiter due to various error sources in the GPS system, 2) frequencies which 
are multiples of once per 12 hours are caused by orbit errors of the GPS satellites 
and 3) orbits for Aristoteles are computed once per day. Therefore, in the worst 
case, frequencies are omitted, or at least down weighted by a certain factor, at 
k + m(u; e /u> 0 ) cpr where \k\ < 5 and m < 5. 

Some remarks 

There are some characteristics of TV -1 due to the choice of the Qyy matrix and the 
structure of the observation equations. In the following it is assumed that the design 
matrix A only consists of observation equations which are computed with the same 
values for a ej r and \i and that Q yy is defined for only one observation type. If one 
assumes white noise then (Tq in (3.23) is a scaling factor for a unit matrix since 
Qyy = ctq/. Accordingly: 

N- 1 = (A T [<TqI)~ 1 A)~ l = (t 0 (A t A)~ 1 (3.24) 

which shows that N~ 1 is simply scaled by parameters determining <Tq [which are 
the instrumental accuracy 07, the sampling time At and the mission duration T). 
Secondly the problem of variation of r in the error analysis is predictable since all 
columns in the A matrix are scaled by a factor (a^/r)^ 1 . As a result any variation 
of these parameters is nothing more than post multiplication of the orginal A matrix 
by a diagonal matrix D containing on the main diagonal values scaling the columns. 

A* = AD => 

- D t (A t Q^A)D => 

(iv*)- 1 = - d-'n-'d- 1 

showing that N _1 is simply pre and post multiplied by D -1 . A diagonal element at 
row i and column i of the inverted normal matrix becomes: 

= NS'D? 
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indicating that a second inversion of N may be avoided. 

The drag problem 

The drag problem in relation to the instrument scaling is much harder to simulate 
in an a priori variance spectrum. Atmospheric density data that is available comes 
from Cactus [elliptic orbit, perigee at 270 km], Dynamics Explorer 2 [elliptic orbit, 
perigee at 270 km] and Atmosphere Explorer- C [circular orbit at about 250 km], cf 
Touboul et al (1990). They show that the atmospheric density strongly depends on 
1) the geomagnetic index K p , 2) the latitude [since fluctuations are a factor of 2 to 
3 larger at high latitudes than at equatorial latitudes] and 3) whether density data 
is taken at night or during the day. An important conclusion is that fluctuations are 
significantly smaller between 0.1 and 0.01 Hz [580. ..58 cpr] than between 0.01 and 
0.005 Hz [58. .27 cpr] and below [< 27 cpr]. 

The easiest way to simulate a drag problem is to apply high pass filtering to the 
gradiometer spectrum above 27 cpr. This is pursued in one of the simulations at 
the end of this Chapter. This simulation is supposed to represent a worst case drag 
situation for Aristoteles as it denies the existence of any lumped coefficient below 
27 cpr whereas it is more likely that a degraded gravity gradient signal is observed 
in this frequency band. 

Constrained least squares solutions 

Some of the results that will be discussed at the end of this Chapter depend on a 
priori covariance information for the unknowns involved in the problem. Consider 
the constrained least squares problem: 

y = Ax + ei (3.25) 

c = Iz + €2 (3.26) 

where eq.(3.26) are constraints to direct the unknowns x to some a priori vector c 
and where I equals to a unit matrix. The solution for this problem is: 

i = (A T Q^A + K- 1 )~ 1 (A T Q^y + c) (3.27) 

where K equals to the a priori covariance matrix of the constraints c in the problem. 
The K matrix describes the a priori covariances of the unknowns x which are sup- 
posed to be centered on the constraints c. Sometimes the constrained least squares 
problem for c = 0 is referred to as a hybrid norm minimization [or collocation] 
problem since (3.27) is the minimum of: 

I T Q-Zz+x T K- l x (3.28) 
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cf (Schrama, 1989a). Here a priori information is used for the potential coefficients 
where K is assumed to be diagonal. The diagonal contains for the a priori r.m.s. 
per coefficient: 

n. = \ji (3.29) 

cf (Kaula, 1966b). In the algorithm, application of a priori information of the un- 
knowns is performed by adding <jj~£ to the diagonal elements of the least squares 
normal matrix. 


3.2.3 Presentation 


After inversion of the normal matrix the so-called variance per coefficient per degree 
6f is computed as: 


' O’ 2 ( C,ra ) + 0-3 ( 5lm ) 

* 1 - 1 ^ 57T1 


(3.30) 


where a 2 [Cimj and a 2 (5/ m j symbolize the diagonal elements of the inverted nor- 
mal matrix at the location of Ci m and 5| m . The following conversions of Si exist in 
order to obtain the degree variance spectra of geoid heights: 


6?(N) = altf(a)6f 


(3.31) 


and gravity anomalies: 


(3.32) 

where /3i(a) is a smoothing operator with a determining the block spacing on a 
sphere for the degree variance expressions of gravity anomalies and geoid heights as 
is described in (Katsambalos,1979). 


3.3 Results 

3.3.1 Gradiometer only results 

Figure 3.1 shows 6i for all gradiometer components up to degree and order 90 in a 
so-called ideal case for Aristoteles. This analysis shows that T uu is consistently the 
most valuable component followed by the off-diagonal terms T uu , T uw and T vwy and 
the remaining diagonal terms T vv and T ww * 

The effect of a limited bandwidth of the gradiometer due to 1) orbit errors, 2) 
attitude problems, and 3) thermal noise effects is shown in figure 3.2. This analysis 
shows that limited bandwidths of the gradiometer seriously affect the outcome of 
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an error analysis. Especially the lower degrees and orders of the gravity solution 
deteriorate rapidly when limiting the lower end of the noise spectrum. 

It was recognized for the gradiometer and the GPS position observation equations 
that the choice of the inclination of the orbited plane plays an important role. This 
problem is illustrated in figure 3.3. It was found that these results strongly depend 
on the value of L [here 120] in the error analysis which is also confirmed by error 
analyses of the Aristoteles gradiometer done by (Koop et al.,1989). 

3.3.2 GPS only results 

Typically the mean r.m.s. per coefficient per degree [6/] derived from gradiometer 
observation equations shows a weak improvement in the lower degree and orders 
[especially in the case where bandwidths are restricted], an optimum sensitivity 
[a minimum] near l & 70 and an exponential deterioration beyond this point as is 
shown in figure 3.2. The best case ”GPS only” results appear to show that the lowest 
degrees and orders are most sensitive followed by a steady exponential deterioration 
toward higher degrees. Examples for GPS on Aristoteles, Gravity probe-B, and 
Topex, are shown in figure 3.4. 

These simulations indicate that gravity fields can be improved up to degree and 
orders around 35, 55 and 85 from GPS derived position information on Topex, GP- 
B and Aristoteles since Si intersects K aula’s rule of thumb at these degrees. The 
corresponding cumulative 1° r.m.s. values for geoid heights and gravity anomalies 
are shown in figures 3.5 and 3.6. The accuracies in terms of geoid heights and gravity 
anomalies look very promising especially in the lower degrees and orders. 

The best and worst cases for GPS on Aristoteles are shown in 3.7 indicating a 
deterioration in the lower degrees and orders; in the worst case a priori standard 
deviations for lumped coefficients at \k\ < 5 and m < 5 are upgraded by a factor 
1000. 


3.3.3 GPS only results, constrained least squares approach 

Figure 3.8 shows Si for Topex, using a least squares approach, assuming / = 65° 
[case 1] and I = 90° [case 2]. The results for I = 65° indicate that a gravity model 
solved from Topex GPS data alone should not exceed / % 15 where it intersects 
Kaula’s rule of thumb. 

In contrast to this result a Topex type of satellite at I — 90° would allow to solve 
for a gravity model up to / « 35. Case 3, 4 and 5 in figure 3.8 show the I — 65° 
results now adding the matrix aK~ x [a is a regularization factor for weighting a 
priori information on the coefficients] to the normal matrix for a = 1, a = 0.1 and 
a = 0.01 respectively. We conclude that: 

1. The results for I — 90° are in any case preferable to those at other inclinations. 

This configuration allows one to solve for a gravity model entirely from one 


observation type coming from one satellite. This is a unique situation since 
most satellite gravity models developed up to now are always "assembled” 
from several orbiters at various inclinations, heights and eccentricities as is 
described for the GEM-T2 model in (Marsh et al.,1989), 

2. The results for I — 65° show that some a-priori information is needed to obtain 
a gravity solution comparable to the I = 90° results since a cross-over occurs 
with "Kaula’s rule of thumb” at / = 35 between a = 0.1 and a = 0.01, 

3. The a = 1, I = 65° collocation solution shows the same characteristics as the 
solutions presented by (Pavlis et al.,1989), Si can not intersect the a priori 
signal curve [a natural property of a Wiener/Kolmogorov type of estimator] 
and it appears that the signal to noise ratio is greater than or equal to 1, 

4. As a direct result of the previous two statements one can conclude that the 
signal to noise ratio above degree 25 is dominated by the choice of a. 

3.3.4 Gradiometer combined with GPS 
Least squares solutions 

The most promising results in terms of the mean r.m.s. per coefficient per degree 
[<fy] are obtained by combining the GPS and gradiometer observation equations as 
is shown in figure 3.9 for / up to 300. In this figure case 1 is the gradiometer 
only solution for Aristoteles, case 2 is the best case GPS solution and case 3 is the 
combination of both solutions. It is estimated that the signal to noise ratio for such 
solutions become equal at degree and order 240. 

Figure 3.10 and 3.11 show the results of gradiometer only and gradiometer with 
GPS solutions converted to point, 1° and 5° mean cumulative geoid and gravity 
anomaly errors. We conclude that: 1) GPS and gradiometer derived solutions are 
complementary, 2) errors in geoid heights are governed particulary by uncertainties 
in the lower degree and orders of the gravity field, 3) the original Aristoteles mission 
objectives [e(Aj) < 5 mgal and e(N) < 10 cm at A = 100 km] are hard to meet 
[or maybe even impossible to meet] without the availability of GPS as a tracking 
facility for Aristoteles and that 4) without the availability of GPS the objectives are 
easier met for A g than N. 

Hybrid norm solutions 

Figure 3.12 shows the results in terms of the mean r.m.s per coefficient per degree 
obtained by combining the GPS and gradiometer observation equations for / up to 
360 pursuing the hybrid norm approach where the full K~ l matrix is added to the 
normal matrix. 
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Curve 1 in figure 3.12 may be considered as a worst case Aristoteles solution at 
I — 96.33° assuming that the atmospheric drag problem prohibits the gradiometer 
to measure any lumped coefficient below 27 cpr. However this solution also depends 
quite heavily on a GPS solution to L=120. Similar solutions to lower L in the GPS 
part revealed unsatisfactory large discontinuities in Si whereas the degree and order 
120 case seemed to be an optimum although still some small jump can be seen. 

Curve 2 in figure 3.12 is a best case Aristoteles solution assuming that below 
27 cpr deteriorated gravity gradients are observed. A simultaneous solution already 
gave satisfactory results when GPS observation equations are added to L = 80 [a 
small jump is still observed at L = 80] and assuming a hybrid norm solution. In 
general one may conclude that this procedure results in a somewhat stronger gravity 
field solution between 1=15 and 120 than the previous case. 

However both solutions show that the inclination problem [a slightly non-polar 
orbit] rad the bandwidth problem may be avoided by adding GPS observation equa- 
tions up to sufficient high degree and order and pursuing a hybrid norm approach. 
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Figure 3.1: Behavior of the gradiometer components T UUI T vv and T ww and T uv , T uw 
and T vw . [h = 200 km, e = 0.001, I = 90°, no bandwidth restrictions, a sampling 
time of 4 s, a mission duration of 6 months, 0.01 E instrument precision, least squares 






Figure 3.2: Effects of limited bandwidths of the 2-axis gradiometer on Aristoteles. 
Case 1: lumped coefficients are considered for 4 < ft < oo, below 27 cpr a 1//3 
behavior is assumed in the observation noise spectrum. Case 2: assuming a flat 
noise spectrum for 4 < /3 < oo. Case 3: no bandwidth limitations. Common 
parameters used in all cases are: h — 200 km, e = 0.001, I = 90° and a mission 
duration of 6 months, 0.01 E instrumental noise, 4 s sampling time, for T uu , T ww 
and ./ tuu • 
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Figure 3.3: The r.m.s per coefficient per degree derived from T uu for I — 90°, 93° 
and 96° at h — 200 km. The mean r.m.s. per coefficient per degree is computed 
using a least squares approach for T uu assuming 4 < /3 < oo, a 1//3 behavior below 
/3 m i n = 27 cpr, a sampling time of 4 s and mission duration of 6 months. 




Figure 3.4: The r.m.s. per coefficient per degree using GPS radial, cross - and 
along track variations for h = 200, 600 and 1300 km at J = 90° assuming 3 cm 
instrumental noise, a sampling time of 1 s and a mission duration of 6 months for 
Aristoteles and mission duration of 24 months for Topex and GP-B. 
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Figure 3.7: The r.m.s. values per coefficient per degree in the GPS best and worst 
case on Aristoteles. In the best case all lumped coefficients are used, in the worst 
case lumped coefficients are down weighted a factor 1000 at \k\ < 5 and |m| < 5. 
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Figure 3.8: The r.m.s. values per coefficient per degree derived from Topex [GPS] 
best case collocation solutions. Here a is a regularization factor for weighting a 
priori information on the potential coefficients. Case 1: I = 65°. Case 2:7 = 90°. 
Case 3: I = 65°, a = 0.01. Case 4: I = 65°, a = 0.1. Case 5: I - 65°, a = 1. 
In all cases we assumed an instrumental precision of 3 cm, a sampling time of 1 s, 
a mission duration of 24 months at an altitude of 1300 km for radial, cross — and 
along track components. 
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Figure 3.9: The r.m.s. values per coefficient per degree for various individual and 
combined GPS and gradiometer solutions. All solutions assume a least squares 
approach and I — 90° and a mission duration of 6 months. Case 1: Gradiometer 
only solution: we assumed 4 < /3 < oo, /3 mtn = 27 cpr, 4 s samphng time and 0.01 E 
instrumental noise for T uu , T uw and T ww . Case 2: GPS only solution: we assumed 
0 < /3 < oo, 1 s sampling time and 3 cm instrumental noise for radial, cross - and 
along track components. Case 3: The combined solution of case 1 and 2. 
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Figure 3.10: The cumulative point, 1° and 5° mean 6i(N) values for combined 
GPS and gradiometer solutions. Dashed: with GPS, solid: without GPS. These 
values are derived from cases 1 and 3 shown in figure 3.9. 
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Figure 3.11: The cumulative point, 1° and 5° mean 6((Ay) values for the combined 
GPS and gradiometer solutions. Dashed: with GPS, solid: without GPS. These 
values are derived from cases 1 and 3 shown in figure 3.9. 
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Figure 3.12: The mean r.m.s. per coefficient per degree for the combined GPS and 
gradiometer solutions pursuing a hybrid norm approach. Curve 1: Gradiometer: 
27 < /3 < oo, L=360; GPS: 0 < /3 < oo, L=120. Curve 2: Gradiometer: 4 < (3 < oo, 
/3 min — 27 (1//3 below L=300; GPS: 0 < /3 < oo, L=80. In both cases we 

assumed a mission duration of 6 months, 4 s sampling time and 6.01 E instrumental 
noise for the gradiometer, while measuring T uu , T uw and T ww . For the GPS receiver 
we assumed a sampling time of 1 s and 3 cm noise in the position estimates. 
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Chapter 4 

Conclusions and 
Recommendations 


The discussion in Chapter 2 has shown that a gradiometer mission at 200 km with 
no drag compensation poses several constraints on the design of the instrument, the 
satellite, the choice of the nominal orbit and the accuracy of attitude restitution. 
In an ideal case, the orbit of a gradiometer satellite should be as low as possible, 
e.g. 160-200 km, near circular and allowing a global coverage of the gravity field 
demanding that I = 90°. However in practice sun- synchronous orbits are chosen 
[I = 96.33° at 200 km] to provide an efficient means of power production using solar 
arrays. 

A 2-axis 0.01 E/\/Hz, gradiometer has been proposed for launch near the end 
of this decade [1996-1998] on a satellite called Aristoteles. The mission objectives 
are to measure the global gravity field in order to obtain geoid heights and gravity 
anomalies to within 10 cm and 5 mgal respectively. 

A treatment of the error sources reveals that orbit errors of several meters appear 
to be no real problem for a 0.01 E/\/Hz gradiometer. The errors caused in the gravity 
gradients are mostly in the low frequency band and can be eliminated by filtering 
the signal below 4 cpr for a 10" 2 E instrument, whereas filtering below 25 cpr is 
needed for a 10 -4 E instrument. 

A static approach to the attitude problem for Aristoteles shows that pitch and 
yaw rotational velocities need to be known to within 5 x 10 -9 and 10“ 8 . rad/s respec- 
tively which poses a severe constraint on the attitude restitution of the instrument, 
cf (ESA, 1989). According to this report modern gyroscopes obtain an accuracy of 
10" 7 rad/s inside the measurement bandwidth which is unfortunately still a factor 
10 to 20 too large for an adequate attitude reconstruction. 

A dynamic approach of the attitude problem shows that rotational velocities, 
and thereby gradiometer signal errors, are decaying at a rate inversely proportional 
to the frequency of the disturbing torque function. A simulation of the dynamic 
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attitude problem discussed in Appendix A shows that the atmospheric torques on 
the non-drag free satellite are probably too high for a 0.01 E gradiometer. 

This means that the static approach already answered the question: highly ac- 
curate gyroscopes possibly combined with star trackers are required for attitude 
restitution. We conclude that additional studies are required to determine whether 
this is feasible by application of modern processing techniques and/or technology 
resulting in an acceptable solution within the budgetary constraints. 

Additionally the problem of scale, coupling and non-linearity errors in the ac- 
celerometers is discussed which allows non- conservative forces such as atmospheric 
drag to degrade the instrument performance _cf Touboul et al. (1990). They con- 
clude that the velocity vector of the spacecraft must be as perpendicular as possible 
to the gradiometer plane. This configuration might require a so-called yaw-steering 
mode of Aristoteles which compensates for a cross-track winds near the poles, see 
also (ESA, 1989). 

In Chapter 3 the results of an analytical error analysis of gravity field parameters 
are discussed assuming various scenarios proposed for Aristoteles [gradiometer and 
GPS receiver], Topex [GPS] and GP-B [GPS]. This analytical technique requires 
a nominal circular orbit having a repeat ratio compatible with the highest degree 
and order of the gravity field. Observation equations for both the GPS and the 
gradiometer part are derived in terms of lumped coefficient equations. The error 
analysis itself is based on variances being the elements of the inverted least squares 
[or hybrid norm] normal matrix which are converted to cumulative mean errors for 
gravity anomalies and geoid heights. 

The error analysis shows that limited bandwidths of the gradiometer of Aris- 
toteles seriously affect the outcome of an error analysis. Especially the lower degree 
and orders of the gravity solution deteriorate rapidly when restricting the lower end 
of the noise spectrum which is related to thermal noise in the gradiometer and e.g. 
atmospheric drag causing disturbing torques on the spacecraft. It was also recog- 
nized for both gradiometer and GPS observation equations that the choice of the 
inclination of the orbited plane plays an important role since the formal errors of 
potential coefficients tend to deteriorate when the inclination is several degrees off 
the polar inclination. 

The most promising solutions for Aristoteles were obtained by combining GPS 
and gradiometer observations. It is shown that 1) GPS and gradiometer derived 
solutions are almost complementary, 2) that errors in recovered geoid heights are 
particulary determined by uncertainties in the lower degree and orders of the gravity 
field, 3) the original Aristoteles mission objectives [e(Asr) < 5 mgal and e(JV) < 10 
cm up to A = 100 km] are hard to meet without the availability of GPS as a 
tracking facility and that 4) without the availability of GPS the objectives are easier 
met for gravity anomalies than geoid heights. A worst case drag simulation using 
gradiometer and GPS observation equations shows that the inclination problem [a 


42 



slightly non- polar orbit] and the bandwidth problem may be avoided by adding GPS 
observation equations to sufficient high degree and order of the geopotential model 
while pursuing a hybrid norm approach. 

The results presented in this study should be interpreted in the sense of an error 
analysis rather than a final solution for the gradiometer/GPS problem. GPS data 
may be processed by numerical techniques as was demonstrated by PavUs et al. 
(1989) for Topex. Gradiometer data could be processed by using the GPS gravity 
field solution simultaneously with the measured tensor components in a least squares 
collocation approach cf (Moritz, 1980) to predict a grid of values on a sphere. An 
actual potential coefficient set, to represent the true nature of the short wavelength 
gravity information, could then be derived by numerical analysis methods using 
orthogonality properties of Legendre functions. This is very similar to the tech- 
nique for deriving gravity field solutions from altimeter data and terrestial gravity 
anomalies cf (Rapp & Cruz, 1986). 
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Appendix A 


Attitude error simulation. 


In this Appendix the dynamic behavior of attitude errors is discussed. The equa- 
tions used in this simulation are (2.20), the non-resonant particular solution of the 
linearized Newton-Euler equations including gravitational torques, and eq. (2.14) 
relating these rotational velocity errors to gravity gradient errors. 

Moments of inertia 

At the time of writing moments of inertia axe not published for Aristoteles. In order 
to avoid the laborious task of computing precise moments of inertia we assumed that 
the principle axes of inertia could be derived from a homogenous cylinder represent- 
ing the satellite’s bus, a plate representing the solar arrays and a thin rod pointing 
forward for the magnetometer boom. Additionally we assumed various dimensions 
and weights of these elements; the total configuration is shown in figure A.l. The 
values found for the moments of inertia are /i = 284.3 [radial anti-Earth pointing], 
J 2 = 476.3 and / 3 = 303.1 kg/m 2 . 

The algorithm 

The algorithm assumes a so-called torque noise level variable [TNL\ which defines 
Ri in eqns. (2.19) as Ri — TNL/I{, This allows to evaluate the derivatives of Q{ 
with respect to time in eq. (2.20) for a given f3 symbolizing the frequency in cpr in 
the torque noise spectrum. The resulting variables 0 X and <9 3 are then substituted 
in (2.14) resulting in sine-cosine expressions for A T u and AT 13 [namely a cos (/3n 0 t) + 
b sin(/3no£)]. The amplitudes c = (a 2 + fr 2 ) 1 / 2 are an indication of the errors in T uu 
and T uw showing the expected 1//3 behavior. Figure A. 2 represents the values of 
TNL on the y-axis and /3 on the x-axis for c = 0.01 and c — 0.0001 E. 

This simulation shows that the lower end of the effective frequency bandwidth of 
T uu [solid line] and T uw [dashed line] is determined by the noise level of the torques 
acting on the satellite. The tensor component T uw shows a minimal frequency ap- 
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Figure A.l: Elements dimensions and weights used for the estimation of moments 
of inertia of Aristoteles. 



Frequency [cpr], revtim=5805.2 sec 

Figure A.2: Maximal allowed torques vs. frequency assuming simulated moments 
of inertia of Aristoteles. 
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Figure A. 3: Symmetry, torques, atmospheric drag 

proximately twice as low as T uu . T ww is not considered in the simulation since it is 
the only tensor component in this configuration that is not modulated by no- [see 
eq. (2.14)] Resonance problems are not plotted in figure A.2, the results are shown 
above 2 cpr. All resonant frequencies occur at 1.08 and 0.29 cpr for Q\ and 0 2 > 
/3 = 1.38 cpr for 0 3 . 

Maximal allowed torque noise level 

In a non-drag free environment the noise level of the torques itself is mainly deter- 
mined by non- conservative forces such as drag acting on the bus and appendages of 
the spacecraft. It is also determined by symmetry of the projected area normal to 
the velocity vector of the spacecraft and the spectral behavior of drag variations. 

An example of symmetry of the projected area, torques and atmospheric drag 
acting on the satellite is displayed in figure A. 3. In this example the torque effect 
due to drag equals to: 

T = 1 pv 2 Cd{hAi - l 2 A 2 ). (A.l) 

For a maximal torque noise level of 3 x 1CT 8 Nm and a 0.01 E gradiometer the 
results displayed in figure A.2 allow the lower end of the bandwidths of T uu and T uw 
to start at respectively 29.58 and 14.80 cpr. At 200 km height p « 3 x 10“ 10 kg jm 3 
and v ~ 7784.3 m/s 2 whereas Cj = 3 for Aristoteles, so that in figure A. 3: 

3 x 10” 8 — ~pv 2 Cd^(l\Ai - l 2 A 2 ) 

requiring that 

A (M* - l 2 A 2 ) < 10~ 6 (A.2) 
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which is a very stringent requirement for symmetry of projected surfaces and their 
distances to the principal axes, even if a so-called yaw steering mode is pursued. 
However Touboul et al. (1990) report that the expected drag fluctuations, having 
wavelengths shorter than 200 s, may contain only 5% of the power of the total 
drag force. Even this assumption might be too stringent for condition (A. 2) since 
it would mean that the uncertainties in /i and I 2 have to remain below 1 to 10 ^m 
which is unlikely taking into account phenomena such as thermal expansion due to 
heating and cooling of the spacecraft. We conclude that attitude restitution must be 
provided by measurements from gyroscopes and star trackers once the gradiometer 
is subject to atmospheric drag. In case a drag free or a shielded gradiometer is 
considered the torque noise level is reduced substantially thereby increasing the 
bandwidth of the gradiometer and relaxing the need for a highly accurate attitude 
reconstruction. 
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Appendix B 


Expressions for gravity 
gradients 


In this Appendix the following coordinate systems and indices are used: 

• U{ = {u, u,u;}: the gradiometer instrument system (u: radial, v: along track 
and w: cross track), 

• r a — {rjUj Q1 u e } and r a = {r,u> 0 , /}: subsets of the total set of orbital parame- 
ters {r,u> 0 ,u; c , /}, 

• x p = {x^y^z}: the geocentric system. 

The following relation exists: 


x = JZ3(— w e ).Ri(— J)Jl3(— w o) 


a 

0 

7 


(B.1) 


where a,/3 and 7 are linearized as a = r + u, f3 = v and 7 = w. Here the poten- 
tial function T is defined in the r a system, see eqns. (3.1) through (3.4), whereas 
derivatives of T are needed in the gradiometer instrument frame it;: 


dT 

dui 


d^dra 
dr a du{ 


(B.2) 


The tensor of second order derivatives in the u t system is obtained by differentiating 
once again with respect to Uj\ 


8 2 T dr a dr b dT d 2 r a 


a 2 r 

duiduj dr a dr b du{ duj ' dr a du{duj 


(B.3) 


53 


To evaluate (B.3) the first and second order derivatives of the r a to the u< system 
are needed. These expressions are derived in the following way: 


1 

d 

\dx p ) 

- 1 dx p 

dui 

d 

8 * 

1 

dui 

d 2 r a 

\ 9xp ] 

- x r 

diiiduj 

.dr a . 

| du{dxij 


d 2 x p dr a drt, \ 
dr a dT\, du{ duj J 


(BA) 

(B.5) 


where the derivatives of x p to r a and of x p to u, are computed from eq.(B.l). 
The formula manipulation program REDUCE, developed by Hearn et al. (1985), is 
used to develop the full tensor of second order derivatives of the form of (B.3). [the 
evaluations of (B.3), (B.4) and (B.5) are rather lengthy] The expressions found which 
are independent of the choice of r a are: [notation: = X r , ^ = T OJ = T e , 

etc.] 


X uu - T rr (B.6) 

T vv = (B.7) 

r* t 

T ww — -T rr -Too T r (B.8) 

r* t 

Tuv = -Tro-^To (B.9) 

v r* 

The terms T uw and T vw depend on the choice of r a : 


and 


T — 

■L UW — 

T uw = 

Tvw = 

Tvw = 


r smw 0 r 

- — 7 — r{(-T e - T re ) + cos/(--T„ + T ro )} 

rcosw 0 smi r r 




t* smu/ 0 
1 




r 2 cos 2 u> 0 sin I 

{ cos uj 0 ( Toe + T 00 cos/) + sinu; 0 (-T e + T 0 cos/)} 


(B.10) 

(B.ll) 

(B.12) 

(B.13) 


Multiplication of (B.10) times sin 2 u; 0 and adding (B.ll) times cos 2 u> 0 results in: 


Xuu; 


r 1 1 _ 1 m . COS / f 1 1 ^-- 

{sin I {—-T re + ^T e } + J^ji-Tro ~ ^ T „ }} cosa; 0 + 

{-T r i — ^Ti} sinwc (B.14) 

r r* 
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In a similar way (B.12) times sin 3 u> 0 added to (B.13) times cos 2 uj 0 results in: 


T — 


-1 


cos/ 


{-j-r—fTeo + -^jToo - -^Tl} COSWo + 
r 2 sini r 2 sini t* 

{ j T t + + ^Tol} sinw 0 

r 2 sini r 2 sinJ r 1 


(B.15) 


While using expressions similar to (B.10) through (B.13) Betti and Sanso (1988) in- 
troduced so-called 77 mp functions which are modifications of the original inclination 
functions. These modified functions are needed to avoid singularities at u> 0 = Ar| in 
eqns. (B.10) to (B.13). Moreover it is desirable to obtain expressions where all time 
dependent effects are contained in the derivatives of T in the r B system [or T* in 
Betti and Sanso’s approach]. 

Here we use expressions (B.14) and (B.15) which merely require to multiply a 
Fourier series by sine and cosine terms thereby avoiding to introduce modified incli- 
nation functions which would require to change the existing algorithm to compute 
inclination functions and their derivatives, cf (Schrama, 1989a). A multiplication of 
a Fourier series once by a sine and once by a cosine term similar to the structure of 
eqns. (B.14) and (B.15) is not very complicated, one can show that: 


equals to 


where 


and 


L L 

X X ( A km COS rpkm + B C km sin V’fcm) COS V’fcm + 

k=—L m=0 



{Ak m COS Ip km + B‘km sin V’fcm) sin V’fcm 

(B.16) 


L+l L 

X X Akm cos + Bkm 

k=-L- 1 m=0 

(B.17) 

Akm 


(B.18) 

Bkm 


(B.19) 


= B l'rn = 0 for 1*1 > L. 

(B.20) 
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